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t— 5 \ Abstract 

(N . 

In this paper we discuss the treatment of discontinuities in Smoothed Particle Hy- 
drodynamics (SPH) simulations. In particular we discuss the difference between in- 
Qh 1 tegral and differential representations of the fluid equations in an SPH context and 

how this relates to the formulation of dissipative terms for the capture of shocks 
' and other discontinuities. 

This has important implications for many problems, in particular related to re- 
cently highlighted problems in treating Kelvin-Helmholtz instabilities across entropy 
gradients in SPH. The specific problems pointed out by Agertz et al. (2007) are 
shown to be related in particular to the (lack of) treatment of contact discontinu- 
| ities in standard SPH formulations which can be cured by the simple application of 

C — ■ an artificial thermal conductivity term. We propose a new formulation of artificial 

thermal conductivity in SPH which minimises dissipation away from discontinuities 
g ; and can therefore be applied quite generally in SPH calculations. 
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1 Introduction 



Smoothed Particle Hydrodynamics (SPH) is a Lagrang ian particle method for 
solving the equations o f fluid dynamics (for reviews, see lMonaghanll2005l ; iPrice 
2004J : lMonagIianlll992l ). 
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Whilst SPH is w idely used in a stroph ysics, geophysics and engineering appli- 



cations, recently lAgertz et al.l (120071 ) have suggested that there are "funda- 
mental differences" between SPH and grid-based codes, particularly relating 
to the simulation of Kelvin-Helmholtz instabilities between two fluids of dif- 
ferent densities. Whilst not phrased in terms of "fundamental differences", 
similar problems have been hinted at previously by other authors regarding 
the treatment of large de nsity gradients in SPH, for example in the context o f 



multi-phase calculations (IRitchie and Thomasl . l200ll ; iMarri and White! . 120031 ) 



The aim of this paper is to resolve these issues in the broader context of 
how discontinuities in SPH are treated, starting from an understanding of 
the difference between integral and differential representations of the fluid 
equations in SPH (in particular the continuity equation) and thus the need for 
"discontinuity-capturing" terms where differential representations are used. 

The paper is structured as follows: Our basic SPH formulation is presented in 
§|2] and the treatment of discontinuities is discussed in ^J3j We propose a new 
formulation for artificial thermal conductivity in SPH in §3] which is both effec- 
tive at resolving contact discontinuities appropriately whilst also minimising 
the dissipation of thermal energy gradients elsewhere. The general discussion 
regarding discontinuties is illustrated on shock tube tests presented in §BZD 
on which the effectiveness of the new artificial thermal conductivity formu- 
lation is also demonstrated. The problems relating to simulating the Kelvin- 
Helmholtz instability in SPH are discussed in §4.21 based on numerical tests. 
We demonstrate that, whilst there are indeed numerical issues with resolving 
KH instabilities across contact discontinuities in standard SPH formulations, 
application of our new term very effectively cures the problem. The results are 
discussed and summarised in §51 



2 Standard variable— h SPH equations 



Whilst the standard derivation of the SPH equations fro m a Lagrangian varia- 
tiona l principle has been pre s ented by many authors fe.g.lNelson and Papaloizoul . 
1994 ; iMonaghan and Price] . 2001 ; Springel and Hernquist . 2002 ; Monaghan 



2002 



Price and Monaghan . 2004bl ) . it is instructive to repeat the de rivation 



here. We begin with the Lagrang ian for a perfect fluid of the form (lEckart 
19601 : iMonaghan and Price! l2001[ ) 



-pv 2 - pu(p, s) 



dV, 



where p, v and u are the fluid density, velocity and thermal energy per unit 
mass respectively, the latter assumed to be a function of density and the 
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entropy s. This integral is discretised as a sum over SPH particles by replacing 
the mass element pdV with the particle mass and the integral by a summation, 
giving 



N 
3=1 



The equations of motion for particle i may then be derived using the Euler- 
Lagrange equations in the form 

d f dL\ _ dL _ 



provided that all of the quantities in the Lagrangian can be expressed as a 
function of the particle co-ordinates r and their time derivatives r = v. In this 
way the exact conservation of momentum, angular momentum and energy is 
guaranteed in the resultant SPH equations because of the symmetry of the 
Lagrangian with respect to translations, rotations and time respectively. 

The momentum is given straightforwardly, from Equation (T5]), by 

tt- = miVi. (4) 



The spatial derivatives in the Lagrangian are found by assuming that the 
entropy is constant and thus that the thermal energy can be expressed as a 
function solely of the fluid density, giving 



dL ^ duj 
— — = — > rrij 



dri 



(5) 



where the derivative of thermal energy with respect to density is provided 
by the first law of thermodynamics at constant entropy, dU= —PdV, where 
V = m/p is the particle volume such that the change in the thermal energy 
per unit mass is given by 

P 

du = —dp. (6) 



The density in SPH is calculated by sum according to 

p i = ^m j W(\v i -r j \,h i ), (7) 
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where in the variable smoothing length formulation hi is in turn assumed to 
be a function of pi, in the form 



h = 7](m/ p) 1 / 3 , 



(8) 



where rj specifies the smoothing length in units of the average particle spacing 
(we use T) — 1.2 throughout this paper). The density summation is thus a non- 
linear equation fo r both p and h which we so lve iteratively using a Newton- 
Raphson method (IPrice and Monaghanl . 120071 ). Taking the time derivative of 
the density sum, we find 



dpi _ 1 

dt ~ a 



(9) 



where Qi is a term relating to the derivative of the kernel with respect to the 
smoothing length. The above is an SPH expression for the continuity equation 
in the form 



dp 
~dl 



-pV • v. 



(10) 



The spatial derivative of density is given by 



dp '1 



giving, via 



and ([3]) the SPH equation of motion in the form 



dvi 
~dt 



The most common method for integrating the energy equation in SPH is to 
evolve the thermal energy, which from (jH]) evolves according to 



dui 



dt Vtip\ 



^rrijivi-Vj) -VWijihi). 



(13) 



Alterna tively the total energy e = \ v 2 + u can be used, which from the Hamil- 
tonian ( Monaghan and Price! . l200l[ ) 



has derivative 



— — = > rrij 
dt 3 



PiV 



(14) 
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A third alternative in the case of an ideal gas is also possible to evolve the 
entropy function S = Pj 1 p"> ', which evolves according to 



dS 


_ 7 - 


1 


( du 


~dt 


P 7 ' 


i 


(dl 




_ 7 - 


1 


( du 




P 7 " 


i 


(di 




diss 

= (no dissipation), (15) 

where the subscript (du/dt)diss indicates the dissipative part of the evolution 
of thermal energy. The latter has the advantage of placing strict controls on 
source s of entropy, since S is purely advected in the absence of dissipative 



terms (ISpringel and Hernquistl . 12002 ) . It should be noted in the context of our 
later discussions that the exact advection of entropy is inherently a differential 
assumption since it relies on the fact that du — Pj p 2 dp = for which there is 
no corresponding integral conservation law. 

The reader unfamiliar with SPH should also note that in SPH evolving the 
thermal energy u differs from an evolution involving the conserved total energy 
e or the entropy A only by the timestepping algorithm. This is quite different 
to the situation in an Eulerian code where there are additional differences due 
to advection terms. In SPH an evolution using either u, e or A will conserve en- 
ergy to timestepping accuracy (assuming the terms associated with smoothing 
length gradients have been properly accounted for as described above). 



3 Discontinuities in SPH 



The treatment of flow discontinuities in numerical hydrodynamics has been the 
subject of a vast body of research over the last 50 years, resulting in the devel- 
opment of a wide range of high accuracy methods for shock capturing schemes 
mainly applicable in the context of grid-based codes. Related to this has been 
an understanding of the assumptions necessary for discontinuous solutions to 
the equations of hydrodynamics (ie. shocks and contact discontinuities) to be 
captured by the numerical solution. Thus for example, the difference between 
a finite volume scheme and a finite difference scheme. Whilst both provide dis- 
cretisations of the fluid equations on an Eulerian mesh, a finite volume scheme 
uses a discretisation of the integral form of the equations, whilst the starting 
point for a finite difference scheme is the discretisation of the fluid equations 
in differential form. The problem with the latter is that the assumption that 
the equations are differentiable immediately excludes the possibility of solu- 
tions which have infinite derivatives (requiring additional mechanisms such as 
dissipative terms in order to capture such solutions), whereas these solutions 
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are not excluded in an integral form. We discuss below how this relates to 
formulation of the SPH equations. 



3.1 Density sum versus density evolution 



Useful insight into the difference between integral and differential formulations 
of the fluid equations in SPH may be gained by considering the difference be- 
tween calculating the density by summation using (j7j) and evolving the density 
as a fluid variable using (J9j) . It is often assumed that these are two equivalent 
ways of calculating the density in an SPH calculation since Equation (Q is 
simply the time derivative of Equation (j7]) and thus that the only distinguish- 
ing factor between the two is the cost associated with calculating the density 
by su mmation separate to the evaluation of SPH forces (see, e.g. iMonaghan . 
1997h . 



In the light of our discussion above, we expect the two to differ because the 
density summation represents an integral formulation of the continuity equa- 
tion (that is, in using the density sum we have nowhere assumed that the 
density is differentiable) whilst in taking the time derivative we have assumed 
that the density is a differentiable quantity. The difference between the two 
can be found by considering the continuity equation written in integral form 
and smoothed over the local volume using the SPH kernel, ie. 



d_£ 
~dt 



+ V • (pV) 



W(\r-r'\,h)dV = 0. 



(16) 



Expanding, we have 
d 



— J p'WdV + J V ■ (p'v')WdV' = 0, 



(17) 



where W = W (\r — r'\, h) . The second term can be expanded further using 
V' ■ [p'v'W] = WV ■ (pV) + pV • V'W, (18) 



giving 





— J p'WdV - J p'V ■ VWdV + J V ■ [p'v'W] dV = 0. (19) 



Making use of the antisymmetry of the kernel gradient (ie. VW = — VW) 
and using Green's theorem to convert the last term from a volume to a surface 
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integral, we have 

J p'WdV + J p'v' ■ VWdV + J [p'v'W] ■ dS = 0. (20) 



d_ 

dt 



Replacing the Eulerian time derivative d/dt with the Lagrangian time deriva- 
tive, ie. d/dt = d/dt — v • V, we find 

j t J p'WdV - J p'{v - V) ■ VWdV + J [p'v'W] ■ dS = 0. (21) 



Finally, we can write the volume integrals as SPH sums by replacing pdV with 
the particle mass and converting the integral to a sum, giving 

j t Y.mjWij = Y, m i(vi - v,) • VW y - / [p'v'W] ■ dS. (22) 

3 j 



The above expression clearly shows that a formulation of the SPH continuity 
equation in integral form would involve not only the time derivative of the 
density sum in the form but also an additional term which appears in the 
above as a surface integral. This term in general vanishes (that is, the kernel 
goes to zero at the limits of the integration volume) except at boundaries, or 
equivalently, flow discontinuities. Thus we expect that SPH formulations which 
evolve the continuity equation in the form (Q will differ from formulations 
which utilise the density sum at such discontinuities, the latter of which should 
require no special treatment. That this is indeed the case is demonstrated 
in numerical tests presented in §H In the context of our variable smoothing 
length SPH formulation, it is important to note that a true "integral form" of 
the density sum is only obtained when the smoothing length is also obtained 
directly from the summation via iteration of the smoothing length-density 
relation discussed in §2] (as opposed to simply evolving h separately using a 
differential form of the continuity equation). 

Regarding discontinuities in other variables (apart from density), the appear- 
ance of surface integrals also provides some insight into where discontinuities, 
e.g. in velocity "go missing" when deriving the SPH equations from a La- 
grangian presented in §2j For example, in using the Euler-Lagrange equations 
([3]) we have implicitly assumed that the variation in the action vanishes at the 
surface of the integration volume (that is, certain surface integrals involving 
the Lagrangian are assumed to vanish). An assumption of differentiability is 
also apparent from the fact that the Euler-Lagrange equations contain deriva- 
tives with respect to particle coordinates and velocity and can only be derived 
in this form by assuming that the variation in the action vanishes at the sur- 
face of the integration volume. The reader will thus note that in our SPH 
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derivation in £J2] a differentiated (and thus assumed different iable) version of 
the density sum was used in finding the equations of motion, leading to the 
inevitable consequence that, whilst the continuity equation can be solved in an 
integral form using the density summation, the SPH momentum and energy 
equations derived above are clearly differential. 



3.2 Artificial dissipation terms 



The discussion above leads to an obvious corollary, namely given that discon- 
tinuities have "gone missing" from the SPH by the assumption of differentia- 
bility, how should they be recovered in the numerical solution? The simplest 
approach is to add dissipation terms to the SPH equations which diffuse dis- 
continuities on the smoothing scale such that they are resolved by the numer- 
ical method (and thus no longer "discont inuous" ) . A gene ral formulation of 
such dissipative terms was presented by iMonaghanl (119971 ) in a comparison 
of SPH to grid-based codes incorporating Riemann solvers. Whilst the usual 
approach taken in S PH is to simply a dd an artificial viscosity term to the 
momentum equation. IMonaghanl (119971 ) noted that, by analogy with Riemann 
solvers, the evolution equation for every conservative variable should contain a 
corresponding dissipation term in it's evolution related to jumps in that vari- 
able, le ading naturally to formulation s of dissipative terms for ultra-relativistic 
shocks ( Chow and Monaghanl. ll997T) and for Magnetohydrodynamics (MHD) 



f lPrice and Monaghanl . l2004al . 120051 ) in SPH. 



3.2.1 Hydrodynamics 



For a non-relativistic gas the dissipation terms for the evolved variables in 
conservative form (namely the conserved mome ntum and energy p er unit mass, 
v and e = \v 2 + u respectively) take the form (IMonaghanl . Il997l ) 




where the bar over the kernel refers to the fact that the kernel must be sym- 
metrised with respect to h, ie. 

VW^ = X - [VWijihi) + VWijihj)} , (25) 
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and the energy variable e* = 2 av sig( v i ' + a uV^ ig Ui refers to an energy 
including only components along the line of sight joining the particles with 
different parameters (a, a u ) specifying the dissipation applied to each compo- 
nent. The choice of signal speed v S i g is discussed below ( §3.2.31) . Note, however, 
that in this paper we have deliberately distinguished between the signal ve- 
locities used for the kinetic energy term v s i g and that used for the thermal 
energy term (v^ ), for r easons that will become clear. This differs from pre- 



vious formulations (e.g. iMonaghanl 119971 ; iPrice and Monaghanl l2004al . 120051 ) 
which have assumed that the same signal velocity is used to treat jumps in all 
variables. 



Equation (|23j) in the IMonaghanl (119971 ) formulation provides an art ificial vis- 



cosity term similar to earlier SPH formulations (e.g. lMonaghanlll992l - the two 



formulations differ only by a factor of h/jjjj]) . Equation (I24p is more inter- 
esting, since (as discussed by Monaghan 19971 ) it shows that the evolution of 
the total energy should contain not only a term relating to jumps in kinetic 
energy (ie. the thermal energy contribution from the viscosity term) but also 
a term relating to jumps in thermal energy. This is more explicitly obvious 
if we consider the evolution of the thermal energy resulting from the above 
formulation, ie. 



du 



de 



dv 
~dt' 



(26) 



which, using ( l23l) and ( l24|) gives 



rrii 



diss j P"*-! 

+ OL u V U sig {Ui - Uj) \ Yij ■ ViWiy 



(27) 



The term involving (ui — Uj) provides an artificial thermal conductivity which 
acts to smooth discontinuities in the thermal energy. The need for such an 
artificial thermal conductivity contribution in order to resolve discontinuities 
in thermal energy is almost universally ignored in SPH formulations. 

The effect of applying different t ypes of dissipation to s pecifi c discontinuities 
is discussed in the M HD case by IPrice and Monaghanl (120051 ) and in the hy- 
drodynamic case by iPricd (120041 ) . The point made in these papers is that 
every physical discontinuity requires an appropriate treatment. For example 
in hydrodynamics, shocks are treated by the application of artificial viscosity 
terms but accurate treatment of contact discontinuities requires the addition 
of artificial thermal conductivity to treat the jump in thermal energy. In the 
MHD case discontinuities in the magnetic field are treated separately by the 
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application of artificial resistivity. We discuss the hydrodynamic case in more 
detail below and in the shock tube tests presented in §4.11 In §4.2l we show how 
these results have a direct bearing on the problems encountered when trying 
to simulate Kelvin-Helmholtz instabilities across density jumps in SPH. 



3.2.2 Interpretation of dissipative terms 



The dissipation terms introduced by IMonaghanl (119971 ) can be interpreted 
more generally as "discontinuity capturing" terms. Interpreted as such, for 
any conservative variable (ie. such that J2j m jdAj/dt = 0) that is evolved 
via a differential equation one would expect to add a dissipation term of the 
general form (for a scalar quantity A) 



'dAj 
dt 



E 



Ifl; 



a A v. 



sig 



diss 



{Ai - Aj)r i3 ■ VWi 



Kl 



(2* 



where is a parameter of order unity specifying the amount of diffusion 
to be added to the evolution of A. The interpretation of ((2SJ) ca n be seen by 



consi dering the SPH expression for the Laplacian in the form (e.g. lBrookshaw 



19851 ) 



{V 2 A\ = 2]Tm J (A Ai) pL (29) 



Pi 



where the scalar function Fij is the dimensionless part of the kernel gradient 
such that VWij = r^Fy and thus • VWy = F^. We then see that (1251) is 
simply an SPH representation of a diffusion term of the form 



— w r?V 2 A (30) 




with a diffusion parameter 77 proportional to the resolution lengthPI 

r] a av sig \rij\. (31) 



3.2.3 Choosing the signal velocity 



In previous formulations (IMonaghanl . 119971 ; iPrice and Monaghanl . l2004bl . l2005l ) 
the signal speed v s i g used in both the artificial viscosity and conductivity 



1 Note that whilst the resolution length appears as the particle spacing, this is 
similar to the smoothing length since within the kernel radius \rij\/h < 2 
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terms is chosen to be an estimate of the magnitude of the maximum signal 
velocity between a particle pair, an estimate for which (for non-relativistic 
hydrodynamics) is given by flMonaghanl . 119971 ) 



'V 




(32) 



where c is the sound speed and generally (3 = 2. However, whilst using a 
signal velocity based on the sound speed and relative particle velocities is 
appropriate at shocks (which travel at the sound speed and involve strong 
compression), it is not clear that the same signal velocity should be used to 
treat contact discontinuities (where there is no compression and the motion is 
at the post-shock velocity). A good example is to consider the simplest case of 
two regions with different densities and temperatures in pressure equilibrium. 
Applying artificial thermal conductivity using a signal velocity proportional to 
the sound speed would result in a steady diffusion of the initial discontinuity 
in thermal energy, which as t — > oo would have completely eliminated the 
temperature gradient between the two regions. 

A much better approach suggested by the shock tube results discussed in §4.11 
is to apply artificial conductivity only in order to eliminate spurious pressure 
gradients across contact discontinuities. In order to do so we require a signal 
velocity which vanishes when the pressure difference between a particle pair 
is zero. We propose the following 



which is constructed to have dimensions of velocity and to be zero once pres- 
sure equilibrium is reached. We find that this is a very effective approach to 
introducing artificial thermal conductivity into SPH in a controlled manner to 
appropriately treat contact discontinuities without the side-effect of unwanted 
diffusion elsewhere. This is particularly the case in the Kelvin-Helmholtz tests 
discussed in §4.21 

3.2.4 Reducing dissipation away from discontinuities 

The key problem with using dissipative terms for capturing discontinuities is 
that such terms also tend to dissipate gradients which are not purely discontin- 
uous. This is a particular problem in relation to artificial thermal conductivity, 
since whilst shocks are continually steepened by the propagating wave, a gra- 
dient in thermal energy, once diffused, will remain diffused forever. The art 
is therefore to come up with well-designed switches that turn the dissipation 
terms off away from discontinuities. 



v 



Pa — P 



(33) 
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In this paper we adopt the artificial viscosity switch suggested by lMorris and Monaghan 



( 119971 ). where the viscosity parameter a is different for every particle and 
evolved according to a simple source and decay equation of the form 



dt 



such that in the absence of sources S, a decays to a value a min over a timescale 
r. The timescale r is calculated according to 

n = (35) 



where h is the particle's smoothing length, v sig is the maximum signal propaga- 
tion speed for particle i (ie. the maximum over pairs involving i of the pairwise 
v S ig defined in Equation [321) and C is a dimensionless parameter which we set 
to C = 0.1 which means that the value of a decays to a m i n over ~ 5 smoothing 
lengths. In general we also impose a maximum value of a max = 1 throughout 
the evolution. In the dissipation terms f[2"3"j) and f[2"T|) the average value on the 
particle pair a = 0.5 (a j + <x,) is used to maintain symmetry. 

The source term S is chosen such that the art ificial dissipation grows as th e 



particle approaches a shock front. We use, as in lMorris and Monaghan! (119971 ) 



S = max(-V • v,0), (36) 



such that the dissipation grows in regions of strong compression. 



A similar switch for the a rtificial th ermal cond uctivity was introduced by 
Price and Monaghan ( 20051 ) (see also Price 2004 ). where the controlling pa- 
rameter a u is evolved according to (1511) with the minimum value a u>min set to 
zero and a source term based on a second derivative of the thermal energy, 



hj\V 2 u\i 
y/Ui + e 



(37) 



where h is the smoothing length, e is a small number to prevent divergences for 
small u and the second derivative term is computed using the standard SPH 
formulation for the Laplacian (Equation 1291). Note that the decay timescale r 
in this case is kept the same as for the viscosity, ie. using ( l32i) . In this paper 
we find that the combination of our new v^ ig and the above switch are very 
effective at turning the conductivity off away from discontinuities (in fact with 
the new v™ ig the switch is almost unnecessary). 
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3.3 Alternative approaches 



Ritchie and ThomasI (bOQlj ) (hereafter RTOl) suggested an alternative ap- 
proach to dealing with problems with contact discontinuities in multiphase 
calculations based on a smoothed estimate of pressure. In their formulation 
the SPH force equation takes the form (modified here slightly to assume an 
adiabatic equation of state and in the form of the symmetry of the kernels 
with respect to the smoothing length) 



alt 



7)E 



11 • ii' 



(38) 



which, translated, is an SPH form of 



alt 



VP P m 
+ -VI 

P P 



(39) 



The mean density (p) used in the force term is derived from a smoothed 
pressure estimate, in the form 

(7 - l)Ui Ui 



The motivation behind this formulation was to compute the pressure force 
without the SPH density explicitly appearing in the equations, in order to 
better handle pressure profiles across strong density gradients. We compare 
the results of this formulation with the standard SPH equations (and to our 
formulation using artificial thermal conductivity) in §4.21 (note that for the 
comparison in this paper we use the smoothing length calculated by iteration 
with the usual SPH density estimate as described in $2j). We indeed find an 
improvement in pressure continuity at discontinuities using their formulation, 
though at the expense of considerable particle noise at the interface. 



Marri and White! (120031 ) have also proposed a modified SPH formulation for 
multiphase flows, in the form of several somewhat ad-hoc criteria for exclud- 
ing particles from one another's neighbour lists. It is however difficult to see 
how their method can be adopted into a consistent Lagrangian formulation 
of the SPH equations, particularly since the total energy and/or momentum 
could arbitrarily change between timesteps as particle pairs are excluded (or 
not) from the calculation. Furthermore the first of their exclusion criteria is 
that the density contrast should be greater than 10, whereas problems with 
resolving the Kelvin-Helmholtz ability occur for much smaller density ratios 
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(as demonstrated in §4.21) where their proposed modifications would have no 
effect. Thus we do not consider their formulation any further here. 



4 Tests 



Whilst there are many problems for which the incorpora tion of an artificial 
therm al conductivity in SPH is a crucial requirement (e.g. iRosswog and Price 
( 120071 ) find significantly improved results on Sedov blast wave tests when such 
a term is incorporated), conciseness limits us to the consideration of two par- 
ticular cases. We first consider a one dimensional shock tube problem in order 
to illustrate the ideas presented in this paper, before applying the method to 
simulations of the Kelvin-Helm holtz instabili t y acro ss a contact discontinuity 
(similar to those performed by Agertz et al.l (120071 ) ). to show that the prob- 
lem highlighted by these authors is very effectively cured by our new artificial 
thermal conductivity formulation. 



4-1 One dimensional shock tube problem 



For the first test we consider a one dimensional shock tube problem where 
an initial discontinuity in pressure and density is set up at the origin. We 
set up the problem with conditions to the left of the discontinuity given by 
(Pl, Pli v l) = (1-0, 1.0, 0.0) and conditions to the right given by (pn, Pr, vr) = 
(0.125,0.1,0.0) with 7 = 5/3. Equal mass particles are used such that the 
particle separation varies according to the density. We choose the spacing to 
the right of the origin as Ar = 0.01 which results in a spacing to the left 
of the origin of — 0.00125, giving a total of 450 particles in the domain 
x = [—0.5,0.5]. Most importantly we set up the problem using unsmoothed 
initial conditions in all variables. Using smoothed initial conditions effectively 
smoothes the initial contact discontinuity throughout the evolution and the 
problems discussed b elow do not become apparent (for a discussion of this 



point, see iPricd 120041 ) . 



It may be noticed in passing that the need for artificial thermal conductivit y 
in this problem has already also been discussed in some detail by lPricd ( 120041 ). 
which we are essentially repeating here because it has a direct bearing on the 
behaviour of Kelvin-Helmholtz instabilities across contact discontinuities and 
thus to the lAgertz et ail (l2007h "problem" . 



To proceed in the order of the discussion given in £j3j we will first begin by 
examining the difference between calculating the density in SPH by direct 
summation in the form and by evolving the continuity equation in the 
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-0.5 0.5 -0.5 0.5 



x x 

Fig. 1. Results of the one dimensional shock tube problem with purely discontinu- 
ous initial conditions using differential forms of the SPH equations. Artificial viscos- 
ity has been applied uniformly and the density (and smoothing length) have been 
evolved using the SPH form of the continuity equation. The contact discontinuity is 
unresolved in the density and results in a spurious 'blip' in the pressure profile and a 
slight overshoot in the thermal energy. The shock is smoothed to several resolution 
lengths by the artificial viscosity term. 

form ([9]). We first consider the evolution in "differential form", that is, using 
the continuity equation. We apply only artificial viscosity in the form (1231) with 
a constant coefficient a = 1 (to keep things simple). The results are shown 
at t — 0.2 in Figure [T], where the SPH particles are indicated by the points 
and may be compared to the exact solution given by the solid line. Whilst 
the shock is smoothed to several resolution lengths by the artificial viscosity 
term, immediately visible is a large 'blip' in the pressure profile around the 
contact discontinuity because of the unresolved gradient in density and a slight 
overshoot in the thermal energy at the same location. 

A significant improvement in the density gradient is obtained by using the 
density summation instead of evolving the continuity equation, shown in Fig- 
ure [2J In the light of the discussion presented in §3.11 we expect no problem 
with the density discontinuities in this case because the density summation 
represents an integral formulation of the continuity equation (put another way, 
the density summation is the exact solution to the SPH continuity equation 
which finds the change in density due to a change in the particle positions). 



15 




-0.5 0.5 -0.5 0.5 



x x 

Fig. 2. As in FigureHJbut where the density has been calculated by direct summation 
which represents an integral formulation of the continuity equation in SPH. The 
density gradient at the contact discontinuity is much closer to the exact solution, 
although there remain problems with the thermal energy gradient and hence the 
pressure. 



That indeed this is the case is demonstrated in Figure [21 which shows that 
using the density summation the density gradient is well resolved across the 
contact discontinuity. However, a spurious blip in the pressure is still apparent 
due to the overshoot in thermal energy. 



The reader who is familiar with so-called "high-resolution shock capturing" 
grid based codes may object at this point that the shock profiles in Fig- 
ure [2] appear excessively smoothed compared to the best grid-based results. 
Although it may be countered that the shock width in a numerical simulation 
is a somewhat meaningless quantity (since the real shock width is many orders 
of magnitud e smaller than the resolution length of an y numerical code), re- 
cent efforts (llnutsukal . 120021 ; ICha and Whitworthl . 120031 ) have shown that such 
methods can also be utili sed in an SPH co ntext. The main objection to doing 
so (raised for example by lMonaghanlll997l ) is that Godunov schemes are more 
difficult to imple ment when the equatio n of state is non-trivial (though meth- 
ods exist - e.g. IColella and Glazl Il985l ). whereas artificial dissipation terms 
are easily applied in all contexts regardless of the degree of complication in 
the physics or the equation of state [for example in an SPH context, the M97 
formulation of the artificial dissipation terms presented above has been readily 
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applied to magnetohyd rodynamics (IPrice and Monaghanl . l2004al . 120051 ) and to 
ultra-relativistic flows (jChow and Monaghanl . 1 19971 )]. 



To show that our key conclusions are not affected by the particular form of the 
"discontinuity capturing terms", Figure [3] shows t he results using the simplest 
(first-order) Godunov-SPH scheme suggested by ICha and Whitworthl (120031 ) 
(that is, where we have simply replaced P, and Pj in the SPH force equation 
with P*, the solution to the Riemann problem treating particles % and j as 
the left and right states, also using P* in the thermal energy equation) in- 
stead of artificial viscosity. The results are almost indistinguishable to those 
using artificial viscosity (Figure [2D apart from a slightly larger smearing of 
the rarefaction (no artificial viscosity is applied to rarefactions in standard 
SPH) but critically also produce a similarly discontinuous pressure at the 
co ntact discontinuity (a simil ar result may be observed in the tests performed 
by ICha and Whitworth 2003 ). This is easily understood since in the simplest 
Godunov-SPH formulation only the jump in pressure has been addressed by 
the scheme and pressure should be constant across a contact discontinuity. It is 
therefore apparent that even in this case som e additional trea tment is required 
to address discontinuities in thermal energy. Ilnutsukal (120021 ) finds similar re- 
sults using Godunov-SPH unless an integral formulation for the momentum 
equation is used based on a polynomial interpolation of density. 



Thus, regardless of the exact form of the "discontinuity capturing terms", 
the key point remains that every evolution equation in SPH resulting from 
a differential formulation requires a term which treats discontinuities in that 
variable. In the present case the missing piece is in the energy equation which 
takes the form of an artificial thermal conductivity. The results of this problem 
using the density summation and including an artificial thermal conductivity 
term of the form (l2"7j) using the new signal velocity (which we expect to diffuse 
the gradient in thermal energy until the pressure is continuous) are shown in 
Figure HI The artificial viscosity has also been applied using the switch given 
in §3.2.41 such that a is a time variable parameter which responds to the 
convergent velocities, demonstrating that the use of such a switch nonetheless 
supplies the necessary dissipation of kinetic energy at a shock front. From 
Figure H] we see that in this case all of the gradients in all of the variables follow 
the exact solution and, most importantly, unlike in all of our previous results, 
the pressure is now continuous (and constant) across the contact discontinuity. 
The continuity of pressure turns out to be of crucial importance in simulating 
Kelvin- Helmholtz instabilities across such discontinuities (see below). 
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Fig. 3. As in Figure [2] (with the density calculated by direct summation) but using a 
Godunov-SPH scheme instead of applying artificial viscosity. The results are almost 
indistinguishable from Figure [2] except for a slightly larger smearing of the rarefac- 
tion in the present case (where artificial viscosity is not normally applied). The 
same problem regarding the thermal energy gradient at the contact discontinuity is 
apparent even in this case. 

4-2 Kelvin- Helmholtz instabilities across a contact discontinuity 



Finally, we consider the problem of Kelvin-Helmholtz instabili ties across a 



densit y jump, on which recent problems have been highlighted by lAgertz et al. 



( 120071 ). We set up the problem in two spatial dimensions, using equal mass 
particles in a periodic box in the domain —0.5 < x < 0.5, —0.5 < y < 0.5. The 
particles are placed on a uniform cubic lattice in separate regions such that 
p = 2 for \y\ < 0.25 and p = 1 elsewhere, where to ensure symmetry in the 
initial conditions we first set up particles in the domain y > and then reflect 
the initial particle distribution across the y = axis. The regions are placed 
in pressure equilibrium with P = 2.5 and we use an ideal gas equation of state 
P = (j—l)pu with 7 = 5/3 such that there is a corresponding jump in thermal 
energy at the contact. As in the one dimensional tests, the initial gradients 
in density and thermal energy are not smoothed in any way (although we do 
first calculate the density by summation before setting the thermal energy to 
ensure that pressure is at least continuous in the initial conditions). Periodic 
boundary conditions are implemented using ghost particles. 
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Fig. 4. As in Figure [2] (with the density calculated by direct summation) but also 
with an artificial thermal conductivity term added using the new pressure-jump 
dependent signal velocity. Artificial viscosity has been applied using the switch 
discussed in §3.2. 41 The density and thermal energy gradients are in this case both 
resolved across the contact discontinuity and most importantly the pressure across 
the contact is continuous, unlike in any of our previous results. 

A cartesian shear flow is setup in the x— direction with velocity v x = 0.5 for 
\y\ < 0.25 and v x = —0.5 elsewhere. Such a configuration is known to be 
unstable to the Kelvin-Helmholtz instability at all wavelengths. In this paper 
we seed the instability at a particular wavelength by applying a small velocity 
perturbation in the y direction given by 

f A sin [-2tt(> + 0.5)/A] \y - 0.25| < 0.025, 

Vy = I (41) 

A sin [2tt(x + 0.5)/A] \y + 0.25| < 0.025. 



where we use A = 1/6 and A = 0.025. 



A textbook anal ysis of the incompressible Kelvin-Helmholtz instability (e.g. 
Choudhuril Il998l ). applicable here since the velocities are smaller than the 



sound speed, shows that the characteristic growth timescale of the instability 
between two shearing layers of densities p and p' is given by 



tkh = 2tt/uj, 



(42) 
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where 
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27T (ppf 



/W 2 I 



(P + P') 



(43) 



For the setup described above, p — 1, p' — 2, v x = 0.5 and = —0.5, giving 
tkh = 0.35 in the units of our simulation for a perturbation of wavelength 
A = 1/6. For the case where the density ratio is 10:1 (ie. p — 1, p' — 10) we 
have tkh = 0.58 for A = 1/6. 

We define the resolution in the simulations according to the particle spacing in 
the least-dense region. For the 2:1 density ratio we show results using particle 
spacing of A = 1/256 and A = 1/512 in this region, resulting in a total of 
97,928 and 393, 160 particles respectively. For the 10:1 case we show results 
using A = 1/256 in the less-dense region, resulting in a total of 359, 604 par- 
ticles. We find that the resolution in the low-density region has an important 
effect on the resolution of the vortex rolls created by the Kelvin-Helmholtz 
instability. Whilst it would be possible to set up the problem using unequal 
mass particles (and thus equal resolution in both regions), it would be rather 
contrived to do so since a density gradient involving equal mass particles is 
the situation which naturally occurs in SPH simulations. 

The results of simulations performed using the setup described and a 2:1 
density ratio are shown in Figure [5] for four different cases (top to bottom), 
with results shown at r KH = 1,2,4,6,8 and 10 (left to right). The four cases 
are, as indicated in the figure: 1) (top) using no artificial viscosity or thermal 
conductivity; 2) using artificial viscosity, applied uniformly with a — 1,(3 — 2; 

3) using the RT01 method with artificial viscosity applied using the switch 
described in §3.2.4t 4) using the usual SPH formulation applying thermal 
conductivity with a u = 1 and our new signal velocity; and finally 5) as in case 

4) but including both thermal conductivity and artificial viscosity (applied 
using switches). Case 2 thus represents the SPH formulation currently most 
widely used (that is, using artificial viscosity with a = 1, /3 = 2) whilst case 5 
represents the formulation we are proposing for general use. 

The differences between the five cases presented in Figure [5] are striking. 
Firstly, in the absence of any dissipation (case 1, top row of Figure [5]), whilst 
there is some growth in the velocity perturbation, mixing between the two 
layers is inhibited by what can only be described as an "artificial surface ten- 
sion" in the dense fluid which results in blobs concentrating into well-defined 
features. This is particularly obvious in an animation of the simulation where 
blobs of dense fluid can be seen to separate and condense into increasingly 
spherical "bubbles" as time advances. Adding artificial viscosity (case 2) does 
not improve the situation, but rather makes things worse by suppressing the 
initial growth in the velocity perturbation for t^h < 6 whilst showing simi- 
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lar surfa ce tension effects at later times. This is in agreement with the results 



found by lAgertz et al.l (120071 ) where reducing the artificial viscosity resulted in 
an increased tendency for layers mix but does not remove the fact that SPH, 
in its standard form, clearly has a problem with Kelvin- Helmholtz instabilities 
across density jumps. 

Use of the RT01 formulation (third row of Figure [5], here shown with viscosity 
applied using the Morris & Monaghan switch) improves things to some ex- 
tent, showing some Kelvin-Helmholtz-instability like features, though at the 
expense of considerable particle noise at the interface (a close-up of which is 
shown in the top left panel of Figure [8]) . 

The results in case 4 (adding artificial thermal conductivity via the method 
discussed in §3.2.41) are substantially different. In this case there is a clear 
development of the classical Kelvin-Helmholtz instability (although still on a 
slower overall timescale than expected from analytic theory), with growth in 
the A = 1/6 mode visible at tkh ~ 2 which is subsequently overtaken by a 
growth in the A = 1/2 mode at t^h ~ 6. The relative growth rate of the 
two modes is, however, in good agreement with theory (that is, a factor of ~3 
difference). Finally, case 5 demonstrated that the dramatic improvement in 
the results found using the artificial thermal conductivity are not significantly 
changed by the turning the artificial viscosity back on (which has been applied 
exactly as in the shock tube test presented in Figure HJ) , largely because the 



Morris and Monaghan! (119971 ) switch ( §3.2.41) is very effective at turning off 



the artificial viscosity where there is no compression. 

An understanding of these results is provided by a plot of the pressure dis- 
tribution in the same set of simulations, shown in Figure [61 In cases 1 and 2 
the boundary between the two fluids is clearly delineated by a ridge of high 
and low pressure. This pressure "blip" at the interface is exactly analogous to 
that observed in the one dimensional shock tube problems (e.g. Figure [2]) and 
is similarly cured by the same solution. Using the RT01 method (case 3) the 
pressure boundary is much less delineated, though considerably noisier, giving 
an improved mixing between the layers compared to the standard SPH formu- 
lation, at the expense of particle noise. Applying the usual SPH method with 
our artificial thermal conductivity term (cases 4 and 5) results in a smooth 
interface and correspondingly good mixing between the layers and a nice KH 
instability. 

To examine the effect of numerical resolution on the results, the five cases 
discussed above are presented at higher resolution (using A = 1/512 in the 
least-dense component) in Figure [3 Comparison of case 1 (top row) shows that 
the artificial surface tension effect is not significantly modified by simply using 
more particles, although the size scale of the "blobs" and "bubbles" of dense 
fluid which break off into the less-dense component are somewhat smaller. 
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Fig. 5. Results of the Kelvin-Helmholtz instability test using a density ratio of 2:1 
and an initial particle spacing of A = 1/256 in the least-dense component. Results 
are shown using (from top to bottom): 1) no artificial viscosity or thermal con- 
ductivity, 2) artificial viscosity applied uniformly, 3) using only artificial thermal 
conductivity applied using our new signal velocity and 4) with both artificial vis- 
cosity and thermal conductivity, where viscosity is applied exactly as in the shock 
tube problem presented in Figure HI 
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Fig. 6. Pressure distribution in the Kelvin-Helmholtz instability test shown in Fig- 
ure Note in particular the regions of high and low pressure delineating the bound- 
ary between the two fluids in the absence of conductivity similar to the pressure 
blip observed around the contact discontinuity in Figure El 
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Fig. 7. Results of the Kelvin-Helmholtz instability test using a density ratio of 2:1 as 
in Figure [5] but here using an initial particle spacing of A = 1/512 in the least-dense 
component. The results are similar to Figure El namely that adding the artificial 
thermal conductivity term gives a dramatic improvement in SPH's ability to resolve 
the Kelvin-Helmholtz instability. 
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Fig. 8. Zoom up of selected panels from Figures [5] and [7] at tkh = 2, highlighting 
the vortex rolls produced using the RT01 formulation (left panels) to a standard 
SPH formulation using our thermal conductivity term (right panels). The RT01 
method works by effectively blurring the discontinuity, which though it helps to 
resolve the Kelvin-Helmholtz instability, also results in considerable particle noise 
at the interface. 

Again, adding artificial viscosity acts to suppress the growth of velocity per- 
turbations (case 2, second row), although in this case some perturbations are 
visible at earlier times compared to Figure El suggesting that the effect of 
viscosity is lessened (which we expect since the artificial viscosity diffusion 
coefficient is linearly proportional to the particle spacing). The RT01 method 
also improves at this resolution, showing clear growth in both the A = 1/6 
mode and the A = 1/2 mode (the latter not well resolved in the lower resolu- 
tion RT01 run in Figure [5]), though the interface remains considerably noisy 
(see close-up shown in the lower left panel of Figure E]). 

In case 4 (applying artificial thermal conductivity using our new v^ ig ) the 
Kelvin-Helmholtz instability is in this case well-resolved for both the A = 1/6 
and the A = 1/2 modes, with well-resolved vortex rolls. The perturbation 
is also much stronger at earlier times (ie. tkh — 1) compared to the lower 
resolution version in Figure [5j Similar results are again found for case 5, in- 
dicating that our proposed "generic formulation" gives good results both for 
this test and on shock-tube problems. The results are also much cleaner than 
those obtained using the RT01 formulation, visible in Figure EJ, though their 
formulation has the merit of being non-dissipative. 

Whilst the tests presented above have been performed, for efficiency, using a 
density ratio of 2:1 between the two fluids, to ensure that our results are not de- 
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Fig. 9. Results of the Kelvin-Helmholtz instability test using a density ratio of 10:1 
with an particle spacing of A = 1/256 i n the least-dense co mponent, shown at the 



(|2007h (t kh = 1/3,2/3 and 



times corresponding to those shown in lAgertz et al. 
1). The SPH calculations are in good agreement with the grid-based calculations 
presented in that paper, although the vortex rolls are slightly less well resolved due 
to the lower resolution employed in the low density fluid. 

pendent on the particular value of the density ratio employed, we have also per- 
forme d a series of calculations using a density ratio of 10:1 (as in lAgertz et al.1 
( 120071 )). The results of this test are shown at early tim es in Figure [9] which may 



be directly compared to the corresponding figure in lAgertz et al. 



(120071 ). In 



particular the results show very good agreement both with the an alytic growth 



timescale and also with the grid-based calculations presented in lAgertz et al. 
( 120071 ). although the vortex rolls are slightly less well resolved in the present 



case due to the slightly lower resolution employed in the low density fluid 
(which differs by a factor of two from the two dimensional equivalent of the 
grid based calculations in the paper because we have used a box twice as large 
in y). The results at later times are shown in Figure[10j Again the development 
of the A = 1/6 mode is clearly seen followed by the A = 1/2 mode, exactly as 
in the 2:1 density ratio case (note that we have scaled the times in units of 
tkh m both cases). Thus we can confirm that o ur results represent a dramatic 
improvement over the SPH results discussed in lAgertz et al.l (120071 ). 



5 Discussion and Conclusions 



In this paper we have discussed several issues related to the treatment of 
discontinuities in SPH which have been shown to ha ve a direct bearing on 
the recent problems highlighted by lAgertz et al.l (120071 ) relating to simulating 
the Kelvin-Helmholtz instability across density jumps. The difference between 
integral and differential forms of the SPH equations was discussed, with the 
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Fig. 10. Results of the Kelvin-Helmholtz instability test, as in Figure [5] but using a 
density ratio of 10:1 applying our new dissipation formulation. 

conclusion that whilst the continuity equation can be formulated in an integral 
form by using the density sum, the SPH momentum and energy equations are 
derived assuming that the equations are differentiable. This leads to the re- 
quirement for diffusion terms in both the momentum and energy evolution to 
capture discontinuities associated with each variable. Whilst discontinuities in 
momentum are treated by the application of artificial viscosity, discontinuities 
in thermal energy (e.g. at a contact discontinuity) were shown to require the 
application of an artificial thermal conductivity term, though this is largely 
ignored in many SPH calculatio ns. This leads to pr oblems at contact discon- 
tinuities such as those found by lAgertz et al.l (120071 ) because of the resultant 
discontinuous pressure profile. 



The reason for the "artificial surface tension" effect was observed to be a 
"ridge" or "blip" in the pressure at the interface between the two fluids in the 
standard SPH formulation (Figure E]). This pressure "b lip" is exactly the ca use 
of the "gap" in the particle distribution discussed by lAgertz et al.l (120071 ). It 
is interes ting therefore to note that in formulations of physical surface ten- 
sion (e.g. iHou et al.lll994j ) the surface tension coefficient is directly related to 
the magnitude of the pressure jump at the interface by the Laplace- Young 
condition, 



[P] = ™, 



(44) 



where [P] = P — P' is the pressure jump normal to the interface, r is the 
surface tension coefficient and k is the interfacial curvature. 

The presence of artificial surface tension in SPH can also be understood as 
a consequence of the exact advection of entropy by SPH particles (which as 
discussed in £j2]is a differential rather than integral conservation law) using the 
following thought experiment (Springel 2005, private communication): Con- 
sider a two-phase distribution of SPH particles where each phase contains 
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particles of different entropy placed in a close box in pressure equilibrium. 
Because entropy is locally conserved there is no difference between a configu- 
ration in which the two phases are separated into distinct spatial regions and 
a configuration in which the two are completely mixed. However if there is 
spurious energy associated with the interface (ie. the pressure is not contin- 
uous) then the unmixed configuration is energetically preferred, leading to a 
tendency for low entropy fluid to form discrete blobs exactly as observed in 
the Kelvin- Helmholtz tests presented in Figures [5] and [71 The effect can be re- 
moved by enforcing pressure continuity across interfaces as we have proposed 
above, which in our case is enforced by adding a small amount of entropy 
mixing (ie. an artificial thermal conductivity) at the interface. Note that this 
argument applies equally to SPH formulations evolving the thermal or total 
energies since such formulations differ only in relation to timestepping pro- 
vided the variable smoothing length terms have been accounted for, and even 
without these terms the errors in the entropy evolution are small. 



To remedy this problem we have presented new formulation for artificial ther- 
mal conductivity in SPH ( §3.2.4|) which acts to equalise pressure at a contact 
discontinuity but which minimises dissipation elsewhere. This term, when ap- 
plied to the Kelvin-Helmholtz problem was found to give excellent results 
in good agreemen t with both analyti c estimates and the grid-based calcula- 
tions presented in lAgertz et al.l (120071 ). We expect this term to be very effec- 
tive at turning off artificial thermal conductivity for hydrodynamic problems 
where pressure gradients always represent non-equilibrium features. However, 
for simulations where other physical effects are included (e.g. gravity) it is pos- 
sible to have pressure gradients that are nonetheless in equilibrium because 
of additional forces in the system. Thus the signal velocity we have proposed 
may require some modification in this case (for example using the "net pres- 
sure difference" accounting for both pressure and gravity between the particle 
pair) although we do not envisage that it would be difficult to do so. 



It is important to note that problems with Kelvin-Helmholtz instabilities oc- 
cur because of entropy gradients which are not treated correctly, not density 
gradients as has often been assumed. Thus a corollary to our results is that 
no special treatment is required in order to capture Kelvin-Helmholtz insta- 
bilities in isothermal flows (as often used in astrophysics) since in this case no 
entropy gradient is present. 



It is worth briefly discussing alternative approaches to the problems described 
above which would be interesting to explore. The first is that, since the prob- 
lem at contact discontinuities appears an "artificial surface tension" effect, a 
fruitful approach might be to try to formulate an "inverse surface tension" 
which counteracts the effect of the surface created by the gradient in particle 
number densit y in the manner of the SPH surface tension formulation recently 
introduced by lHu and Adamsl (120061 ). The advantage would be that this would 
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be a truly dissipation-less approach to tackling contact discontinuities, going 
some way towards actually accounting for the surface integral terms which 
have "gone missing" from the standard SPH formulation by the assumption 
of differentiability. A second approach would be to try to actually calculate 
the surface terms associated with discontinuities in momentum and thermal 



energy (for example iKatd (120011 ) presents a Lagrangian variational principle 
which incorporates discontinuities by retaining sur face integral term s) . Thirdly 
our expectation based on the results presented in llnutsukal (120021 ) is that the 
formulation of Godunov-SPH which accounts for the surface integral terms 
should also give good results for the Kelvin-Helmholtz instability problem 
because of the continuity of pressure across the contact discontinuity. 



In summary, based on the fact that the Kelvin-Helmholtz instability is strongly 
inhibited across density jumps in "standard" SPH formulations (adopting only 
an artificial viscosity term) and the effectiveness with which our proposed 
signal velocity turns off artificial thermal conductivity where it is not needed 
(at least for non-self-gravitating calculations), we suggest that the "generic 
formulation" of dissipation terms presented in this paper should be generally 
adopted as standard procedure in SPH calculations. 



Acknowledgments 



DJP is supported by a UK PPARC/STFC postdoctoral research fellowship. 
I would like to thank the many people who have encouraged me to write up 
the ideas presented in this paper, including Matthew Bate, Giuseppe Lodato, 
Jim Pringle, Thorsten Naab, Markus Wetzstein and Walter Dehnen. Extra 
thanks go to Joe Monaghan from whom I have learnt most of the ideas dis- 
cussed in this paper and to Seung-Hoon Cha for useful discussions regard- 
ing Godunov-SPH. Thanks also to Mike Norman and Volker Springel for 
helpful di scussions. F igures and exact solutions have been produced using 
SPLASH (jPricel . 120071 ). a visualisation tool for SPH that is publicly available 
at http : / / www. astro . ex. ac . uk/p eople / dprice 



References 

Agertz, O., Moore, B., Stadel, J., Potter, D., Miniati, F., Read, J., Mayer, L., 
Gawryszczak, A., Kravtsov, A., Nordlund, A., Pearce, F., Quilis, V., Rudd, 
D., Springel, V., Stone, J., Tasker, E., Teyssier, R., Wadsley, J., Walder, 
R., Sep. 2007. Fundamental differences between SPH and grid methods. 
MNRAS 380, 963-978. 



29 



Brookshaw, L., 1985. A method of calculating radiative heat diffusion in par- 
ticle simulations. Proceedings of the Astronomical Society of Australia 6, 
207-210. 

Cha, S.-H., Whitworth, A. P., Mar. 2003. Implementations and tests of 

Godunov-type particle hydrodynamics. MNRAS 340, 73-90. 
Choudhuri, A. R., Dec. 1998. The physics of fluids and plasmas : an 

introduction for astrophysicists /. New York : Cambridge University 

Press. QB466.F58 C46 1998. 
Chow, E., Monaghan, J. J., 1997. Ultrarelativistic SPH. J. Comp. Phys. 134, 

296-305. 

Colella, P., Glaz, H. M., Jun. 1985. Efficient solution algorithms for the Rie- 

mann problem for real gases. J. Comp. Phys. 59, 264-289. 
Eckart, C, 1960. Physics of Fluids 3, 421. 

Hou, T. Y., Lowengrub, J. S., Shelley, M. J., Oct. 1994. Removing the Stiffness 
from Interfacial Flows with Surface Tension. J. Comp. Phys. 114, 312-338. 

Hu, X. Y., Adams, N. A., Apr. 2006. A multi-phase SPH method for macro- 
scopic and mesoscopic flows. J. Comp. Phys. 213, 844-861. 

Inutsuka, S., Jun. 2002. Reformulation of Smoothed Particle Hydrodynamics 
with Riemann Solver. J. Comp. Phys. 179, 238-267. 

Kats, A. V., May 2001. Variational principle and canonical variables in hy- 
drodynamics with discontinuities. Physica D, 459-474. 

Marri, S., White, S. D. M., Oct. 2003. Smoothed particle hydrodynamics for 
galaxy-formation simulations: improved treatments of multiphase gas, of 
star formation and of supernovae feedback. MNRAS 345, 561-574. 

Monaghan, J. J., 1992. Smoothed particle hydrodynamics. Ann. Rev. Astron. 
Astrophys. 30, 543-574. 

Monaghan, J. J., Sep. 1997. SPH and Riemann Solvers. J. Comp. Phys. 136, 
298-307. 

Monaghan, J. J., Sep. 2002. SPH compressible turbulence. MNRAS 335, 843- 
852. 

Monaghan, J. J., 2005. Smoothed particle hydrodynamics. Rep. Prog. Phys. 
68 (8), 1703-1759. 

Monaghan, J. J., Price, D. J., Dec. 2001. Variational principles for relativistic 

smoothed particle hydrodynamics. MNRAS 328, 381-392. 
Morris, J. P., Monaghan, J. J., 1997. A switch to reduce SPH viscosity. J. 

Comp. Phys. 136, 41-50. 
Nelson, R. P., Papaloizou, J. C. B., Sep. 1994. Variable Smoothing Lengths 

and Energy Conservation in Smoothed Particle Hydrodynamics. MNRAS 

270, 1. 

Price, D. J., 2004. Magnetic fields in Astrophysics. Ph.D. thesis, University of 
Cambridge, Cambridge, UK. astro-ph/0507472. 

Price, D. J., Oct. 2007. SPLASH: An Interactive Visualisation Tool for 
Smoothed Particle Hydrodynamics Simulations. Publications of the Astro- 
nomical Society of Australia 24, 159-173. 

Price, D. J., Monaghan, J. J., 2004a. Smoothed Particle Magnet ohydrody- 



30 



namics I. Algorithms and tests in one dimension. MNRAS 348, 123-138. 
Price, D. J., Monaghan, J. J., 2004b. Smoothed Particle Magnetohydrodynam- 
ics II. Variational principles and variable smoothing length terms. MNRAS 
348, 139. 

Price, D. J., Monaghan, J. J., 2005. Smoothed Particle Magnetohydrodynam- 
ics III. Multidimensional tests and the V • B = constraint. MNRAS 364, 
384-406. 

Price, D. J., Monaghan, J. J., Feb. 2007. An energy- conserving formalism for 
adaptive gravitational force softening in smoothed particle hydrodynamics 
and N-body codes. MNRAS 374, 1347-1358. 

Ritchie, B. W., Thomas, P. A., May 2001. Multiphase smoothed-particle hy- 
drodynamics. MNRAS 323, 743-756. 

Rosswog, S., Price, D., Aug. 2007. MAGMA: a three-dimensional, Lagrangian 
magnetohydrodynamics code for merger applications. MNRAS 379, 915- 
931. 

Springel, V., Hernquist, L., Jul. 2002. Cosmological smoothed particle hydro- 
dynamics simulations: the entropy equation. MNRAS 333, 649-664. 



31 



